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1.  INTRODUCTION 


The  theory  of  isosiasy  has  a  well-documented  history  of  more  than  a  hundred 
years.  Its  fundamentals  are  even  considerably  older  and  can  be  traced  back 
to  Archimedes.  From  geodetic  evidences  some  kind  of  isostatic  equilibrium  had 
to  be  postulated.  Although  in  some  limited  areas  the  Pratt/Hayford  system 
seemed  to  prevail,  the  Airy/Heiskanen  system  is  now  generally  believed  to 
model  the  complex  reality  better. 

The  geodetic  interest  in  isostasy  was  considerably  fueled  both  by  the  access 
to  high  speed  computers  and  the  release  of  a  worldwide  digital  terrain  model. 
Harmonic  coefficients  of  the  topographic-isostatic  potential  have  been 
calculated  by  several  colleagues  up  to  degree  and  order  36,  based  on  5*  «  5* 
mean  topographic  information,  both  in  linear  approximation  (Khan,  1973)  and  in 
the  non-linear  mode  (Lachapelle,  1975).  At  that  time  much  higher  resolutions 
could  not  be  obtained  because  of  missing  data,  but  in  particular  because  of 
excessive  computer  time  requirements.  With  the  design  of  the  fast  Fourier 
algorithm  on  the  sphere  by  Colombo  (1981)  and  the  public  release  of  his 
programs,  high  resolution  models  became  suddenly  feasible  and  have  been 
computed  complete  up  to  degree  180  by  Rapp  (1982a)  based  on  the  worldwide 
digital  topographic  model  of  64800  1*  «  1*  mean  elevations  provided  by  DMAAC 
in  1979  and  using  the  Airy/Heiskanen  model  in  linear  approximation.  Rapp 
claimed  that  correlation  studies  between  his  geopotential  solution  of  1981  and 
the  topographic-isostatic  potential  suggest  a  considerably  deeper  level  of 
compensation  of  about  SO  km  rather  than  the  usually  used  value  of  30  km. 
Rapp’s  linear  approximation  has  been  extended  recently,  taking  into  account 
also  second  suid  third  order  terms,  by  Rummel  (private  communication).  For  a 
compensation  depth  of  50  km  the  agreement  between  the  topography-isostasy 
and  the  geopotential  spectrum  was  very  good  in  the  frequency  range  between 
degree  50  and  150;  however,  rather  poor  results  have  been  obtedned  for  the 
range  between  degree  20  and  50  and,  what  is  more  astonishing,  also  in  the 
high  frequency  range  above  150. 

In  1983  an  "exact"  solution  has  been  produced  by  former  colleagues  at  the 
Technical  University  in  Graz  using  the  setme  DMAAC  data  set,  but 
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unfortunately  a  less  exact  computer  program  which  yielded  >  softly  speaking, 
not  utterly  reliable  results.  Since  this  data  set  is  often  identified  with  my 
name  for  unidentified  reasons,  I  feel  obliged  to  clarify:  it  is  not  my 
product. 

Tscherning  (1984)  compared  existing  solutions,  found  that  the  above 
referenced  data  set  must  be  in  doubt,  criticized  Rapp’s  1982  solution  because 
of  the  too  large  and  therefore,  hard  to  sell  depth  of  compensation,  and 
suggested  the  use  of  much  smaller  compensation  depths  in  the  range  between 
15  and  35  km  which  he  found  "opUmal"  in  terms  of  maximally  reducing  the 
power  of  the  observed  geopotential. 

The  question  of  the  meat  between  the  sandwich  triggered  my  interest  and  was 
a  challenge  to  investigate  that  problem  from  scratch.  The  goal  of  this  study 
was  to  estimate  the  parameter(s)  of  the  most  likely  compensation  model  which 
is  in  best  possible  agreement  with  the  observed  geopotential  and  also 
geophysically  acceptable,  and  the  calculation  of  the  harmonic  coefficients  of 
the  topographic-iaostatic  potential  complete  up  to  degree  and  order  180,  based 
on  that  model. 

Now  you  may  continue  with  Chapter  6,  the  other  part  of  the  sandwich,  and  if 
you  are  turned  on,  you  may  even  find  some  meat  between. 


2.  HARMONIC  ANALYSIS  OP  THE  TOPOGRAPHIC-ISOSTATIC  POTENTIAL 


The  topotfraphic-ieostatic  potential  of  Airy/Heiskanen  type  Tn  ia  defined  as 
the  potential  of  all  mass  disturbances  relative  to  an  ideal  crustal  layer  of 
uniform  density  po  and  thickness  Dt  superimposed  upon  underlying  material  of 
equally  uniform  density,  pi.  Denoting  the  mass  disturbance  by  6p,  the 
potential  Tti  is  well-known  to  be  given  by 

Tti(P)  =  G  III  (P.Q)  MO)  <lv(Q)  (2-1) 

V 

with  G  ...  Newton's  gravitational  constant, 

i(P,Q)  ...  apace  distance  between  P  and  Q, 
dv  ...  volume  element. 

Tti  certainly  harmonic  outside  the  earth's  surface  and  Tu’s  spherical 
harmonic  series  is  certainly  convergent  outside  a  sphere  completely  enclosing 
the  earth.  Outside  that  sphere  we  may  use  the  convergent  series 
representation  of  I”*, 

<-‘(P.0)  =  ZZ  4lr  Pn  (cos  ypg)  (2-2) 

n=o  »  p 

with  r  ...  modulus  of  the  radius  vector, 

fpa  ...  spherical  distance  between  P  and  Q, 

P„  ...  Legendre  polynomial  of  degree  n 

(Heiskanen  k  Moritz,  1967  (HM),  p.  33).  Decomposing  P„  in  terms  of 

n 

P„(cos  Ypg)  =  ^  ^  R„.(P)  R„«(Q)  (2-3) 

with  the  fully  noraalized  spherical  hanonics  R„mi 

yr ^  f  COB  ■  Xp  for  ■  <  0 

2‘"®-(2n+l){^^  P„.(coe  9p)  (2-4) 

•  I  sin  ■  Xp  for  ■  >  0 

...  Kronecker  symbol,  P„^  ...  Legendre  function,  9,  X  ...  polar  distance, 
longitude.),  the  topographic-isostatic  i>otential  is  represented  by  the  harmonic 


senes 


Tti(P)  =  G^V("+‘)  n|fip(Q)r"(Q)R„,(Q)dv(Q)  (2-5) 

m—  ft  ^ 

With  dv(Q)  =  r*(Q)  dr(Q)  d<r(Q) 
and  the  spherical  surface  element  de. 


Let  us  now  investigate  the  volume  integral,  confining  ourselves  to  the  usual 
spherical  approximation.  Then  its  contribution  duo  to  the  topographical 
masses  outside  the  geoid  and  ocean  water  inside  the  geoid  is 

R+H 

//  |tfp(Q)r"+»(Q)dr(Q)H„,(Q)<W(Q),  (2-6a) 

a  r=l« 

and  ita  contnbution  due  to  '.he  isostatically  compensating  masses  is 


R-D 

1/  I  MQ)r''+>(Q)dr(a)R„,(Q)de(a), 

9  r=R— 0+kH 


with  R  .  .  . 
D  .  .  . 
H  .  .  . 

h  •  .  * 


mean  earth  radius, 
depth  of  compensation  level, 
topographical  height  {+), 
ocean  bottom  height  (-), 
compensation  "factor”. 


(2-6b) 


(The  compensation  "factor"  k,  which  is  in  the  non-planar  case  actually  a 
function  sUghtly  dependent  on  H,  will  be  discussed  later.)  In  (2-6a)  the 
function 


=  Po  =  constant  for  topographic  aasses, 
^P  =  Pa  ~  Po  =  constant  for  ocean  aasses. 


in  (2-6b) 

^P  =  Po  ~  Pi  =  constant  for  coapensating  aasses 
is  used  -  density  of  ocean  water).  The  integration  of  (2-6a,b)  with 


4 


4 


respect  to  r  is  straightforward  and  yields 


-  1  R„.(Q)de(Q) 


(2-6a)’ 


(I  -  ir  ^  ^  r  - 


(2-6b)’ 


respectively. 


The  compensation  "factor"  k 


This  concept  of  iaostasy  is  based  on  the  principle  of  mass  balance:  surplus  + 
deficit  z  0.  In  the  simple  Airy/Heiskanen  isostatic  model  this  enables  us  to 
derive  the  geometry  of  the  compensation  surface  (roots/antiroots)  such  that  an 
exact  1  :  1  relation  is  postulated  between  a  topographic  height  H  and  the 
corresponding  compensation  height  kH.  This  mass  balance  principle)  applied  to 
rock  topography  and  its  compensation)  yields 

PoR'[(l  >  D’  -  l]  (Po  -  Pi)(R  -  D)>[l  -  (l  +  =  0)  (2-7) 

and  after  elemenUu-y  manipulations 


kH 

R  -  D 


.  1  -  (i  -  D"’  k.  [(i .  D’  - 1]}’  - 


.  rock;  (2-7a) 


for  the  ocean  and  its  compensation  we  obtain  a  similar  expression) 


R  -  D 


ocean  (2-7b) 


(Pl  -  Po) 


It  is  instructive  to  investij^ate  the  planar  case.  For  this  purpose  we  evaluate 
kH  in  a  power  series  with  respect  to  H, 


and  if  we  let  the  radius  R  become  infinite,  we  obtain 

k  =  -koCH  (2-8) 

with  c^  =  1  for  H  >  0  and  =  1  -  ^  for  H  <  0. 


It  is  obvious  that  the  compensation  factor  k  ia  constant  in  the  planar  model 
and  independent  of  both  H  and  D.  For  standard  densities  po  -  2.67,  pi  -  3.27, 
Pu,  =  1.03  gcm~’,  we  obtain  the  familiar  values  of  -4.45  and  -2.73  for  rock  and 
ocean  compensation,  respectively  (HM  p.l36).  In  the  spherical  model  the 
compensation  "factor"  k  ia  no  longer  a  constant;  it  is  a  function  which 
depends  to  a  small  extent  on  both  the  height  H  and  the  depth  of  the 
compensation  level  D  according  to  (2-7a,b)  with  (2-8)  as  the  leading  term.  In 
order  to  simplify,  we  introduce  "normalized"  heights  H  and  H'  according  to 


(2-9) 


and  the  earth’s  mass  M  by 


M  =  ^  PS  (2-10) 

with  the  mean  earth  density  ~  5.517  gcm~*.  With  this  notation,  the 
topographic-isostatic  potential  in  its  haznnonic  series  representation  is  given 
by 


(2-11) 


with  the  harmonic  (Fourier)  coefficients  of  the  topographic-isostatic  potential 


3.  APPROXIMATIONS  FOR  THE  HARMONIC  COEFFICIENTS  OP  THE 


TOPOGRAPHIC-ISOSTATIC  POTENTIAL 


3.1  Linear  APDroximation 


Given  a  worldwide  digital  terrain  model  (DTM)  representing  the  geooietry  of 
the  solid  earth  surface,  H(Q),  standard  numerical  integration  techniques  have 
been  used  in  the  past  (Khan,  1973;  Lachapelle,  1975).  Due  to  the  excessive 
computation  requirements,  only  low  degree  coefficients  (up  to  degree  and 
order  36)  using  a  5*  x  5*  DTM,  could  be  determined.  A  couple  of  years  ago, 
Colombo  developed  a  fast  Fourier  algorithm  for  the  sphere  which  can  be 
efficiently  applied  if  the  integrand  is  independent  of  any  degree  n. 
Unfortunately,  this  is  not  the  case  with  (2-12b,c).  However,  due  to  the 
smallness  of  H  and  H*  compared  to  1,  a  linear  approximation  seemed  to  be 
justified,  yielding  integrands 


^  [(1  +  H)"-»->  -  1]  =  H, 

^4-3  C<i  +  - 1]  -  »'• 


(3-1) 


which  are  obviously  independent  of  the  degree  n.  Therefore,  in  the  linear 
approximation,  the  FPT  algorithm  applied  to  H  and  H',  respectively,  can  be  and 
has  been  used  (Rapp,  1982),  providing  the  Fourier  coefficients 


^  JJ  c„H(Q)R„«(Q)d^(Q). 

'  (3-2) 

^  JJ  H’(Q)H„.(Q)dv(Q). 

a 

in  a  very  efficient  way. 


In  this  linear  approxioiation  case,  the  product  ChH  is  usually  called  "equivalent 
rock  topography"  which,  loosely  speaking,  trades  in  rock  for  water.  The 
products  />oChH  and  poko'H*  can  consequently  be  interpretated  as  surface 
layer  densities  at  zero  level  and  at  the  level  of  compensation.  Therefore,  in 
the  linear  approxisiation,  the  topographic-iaoetatic  potential  is  repreaented  in 
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terms  of  a  double  layer  potential. 


The  linear  approximation  of  (2-12b)  yields 

TSIiT  5;  a  o«H(a)R„.(Q)<i.(a); 


(3-3a) 


the  linecu*  term  of  (2-7b)  is 


H’  i  -(1  -  D”’koCHH 


leading  to 


-  -  5;  (sfe  -  il"  s:  n  <=«H«»s~(«)<i««»- 


(3-3b) 


Adding  and  according  to  (2~12a)  the  harmonic  series  of  the 

topographic-iaostatic  potential's  linear  approximation  is  given  by 


T„(P)  ^  4=^  tr  ihV  [*  -  (1  -  i)”] 

n  -  - 

•  ^  Km(P)  I;  jJ  CMH(Q)R„.(Q)d«r«l) 


(3-4) 


Exactly  the  same  expression  can  be  obtained  if  rock/ocean  and  compensation 
masses  are  considered  as  single  layers  at  zero  level  and  at  the  level  of 
compensation,  respectively. 

Equation  (3~4)  is  well  suited  to  discuss  the  two  extreme  cases  of  D  =  0  and 
D  =  R  .  It  is  obvious  that  for  D  =  0  the  topographic-isostatic  potential  in  its 
linear  approximation  vanishes  identically.  For  the  second  case,  D  =  R,  the 
topographic-isostatic  potential  reduces  to  the  topographic  potential  with 
vanishing  zero  degree  harmonic  (due  to  the  mass  balance).  Neither  of  these 
two  cases  is  realistic,  as  we  know  from  geodetic  and  geophysical  evidence.  If 
we  consider  "reasonable"  values  of  0  <  D  «  R,  we  may  (at  lecmt  for  low 
degrees)  approxiaiate 


(3-5) 


'  o. .  • . " 


we  observe  that  the  main  terms  in  (3-4)  cancel  each  other  and  obtain  the 
approximation 

Tt,(P)  i  4=^ LI  (n)"  ^  WW  S  f!  c,H(Q)B,.(0)d.(Q) 

(3-6) 


which  can  be  further  simplified  with 

n  ^  1 
2n  +  1  2 

and  the  harmonic  representation  of  the  equivalent  rock  topography, 


(3-7) 


ChH(P)  =  Y1  IZ^neCP)  fc  CMfl(Q)i«,(Q)ds(Q). 

fi  m 


we  find  at  zero  level 


(3-8) 


Tti(P)  =  27iGDpoCmH(P) 


(3-9) 


(cf.  also  HM,  p.  149).  This  approximative  formula  is  a  very  useful  rule  of 
thumb  which  provides  a  fairly  good  estimate  for  a  standard  Airy/Heiskanen 
isostatic  model.  Prom  this  simple  expression  we  conclude  that  the 
topographic-isostatic  potential  is  approximately  linearly  dependent  on  the 
depth  of  compensation  and  also  linearly  on  the  height  of  the  equivalent  rock 
topography. 


3.2  Higher  order  aporoximutiona 


Rapp's  (1982)  argument  that  the  efficient  FFT  algorithm  is  only  applicable  in 
the  linear  approximation,  which  we  have  shown  is  the  double  layer  formulation, 
is  no  longer  valid  if  we  use  the  binomial  expansion  of  the  left-hand  side  of 
(3-1), 


(1  +  H)"+»  -  1  =  ^  (“  j  H’’ 


(3-10) 


(analogous  for  H’),  yielding  exact  expressions  (within  spherical  approximation) 
of  the  harmonic  coefficients  of  the  rock/ocean  topography  and  its  isostatic 


compensation: 


(2n+l)(n+3) 


StI  1  4Tr  JJ  c„H^(0)H„,(Q)de(Q) 

J~‘  a 


(3-lla) 


(2n+l)(iH-3) 


j  )  1;^  jj  CMfl--'(Q)H„.(0)d«r(0) 
J“*  a 


(3-llb) 


There  is  no  reason  why  PFT  could  not  be  applied  to  any  power  of  H  and  H’, 
resi>ectively.  (Note  that  (3-1  la, b)  are  exact.)  But  do  we  gain  anything  by 
employing  FFT,  considering  the  fact  that  2(n  4-  3)  PFT*s  have  to  be  performed 
(one  PFT  for  each  power  of  H  6tnd  S'),  competred  to  twice  a  standard  numerical 
integration?  In  the  case  of  our  earth  we  know  that  H,  H'<<  1;  therefore, 
(3-10)  converges  very  quickly  and  we  can  safely  terminate  the  binomial 
expansion  at  a  very  low  degree  without  committing  a  significant  error  (an 
upper  limit  of  power  j  =  5  is  perfectly  sufficient).  Consequently,  only  a  small 
number  of  PFT'a  must  be  calculated;  we  typically  save  a  tremendous  amount  of 
computer  time  without  sacrificing  accuracy,  and  can  arbitrarily  improve  the 
result  by  adding  another  power  of  H  and  H*  if  we  so  desire. 


3.3  The  second  order  effe 


By  the  second  order  effect  we  understand  the  contribution  of  the  second 
power  of  H  and  H'  to  the  topographic-isostatic  potential.  Because  of  the 
rapid  convergence  of  (3-10),  this  second  order  effect  represents  with  very 
good  approximation  the  difference  between  the  actual  topographic-isostatic 
potential  and  its  double  layer  approximation.  Is  this  difference  sufficiently 
small  to  be  safely  neglected? 


Actually,  there  are  two  kinds  of  second  order  effects:  one  is  due  to  the 
second  power  of  H  and  H'  in  the  binomial  expeuision  (3-10),  the  other  is  due 
to  the  second  power  of  H  in  the  compensation  function  k  =  k(H)  of  (2-7a,b). 


The  second  order  effect  due  to  rock/ocean  topography  is  straightforward: 
using  the  binomial  expansion  (3-10)  in  (2-12b)  and  observing  the  binomi^ 


coefficient 


we  obtain 

‘  (IsTif  h  H'"®’  ■  (S-IZ*) 

For  the  compensation  part  we  develop  (2-7atb)  in  a  series  up  to  the  second 
power  of  Ht 


(l  -  |)*H’  =  -k*CMH{l  +  H)  -  kgcjli*(l  -  §)"*  +  0(i»). 
and  using  again  the  binomial  expansion  (3-10)*  we  get 

-  i)’“'  =  (f  *  i]mcSS*(i  -  2)'*  *  o(i*) 


and  obtain  as  second  order  contribution  of  the  isostatic  compensation 


jS(c) 


3po  1 
PH  (2n+l) 


(*  -  a)"  j;  1/  '«*’[!  *'•=«('  -  i)”  - 

9 

(3-12b) 


Synthesizing  and  according  to  (2-11),  we  obtain  the  sec<»d  order 

effect  of  the  topographic-isostatic  potential. 


•  n 


(3-13) 


with  the  harmonic  coefficients 
fiT„.  :=  «?<;*)  + 

=■  ^  (di)  h  fJ'-S'di  *  i)*(»  -  Si’ll  ‘.'"(i  -  SI”  -  *1)""-®'- 

(3-14) 

Using  the  approximations  (3-5)  and  (3-7),  this  odd-looking  expression  can  be 


considerably  simplified  if  we  admit  an  error  on  the  order  of  20X,  which  can  be 
accepted  since  the  second  order  effect  is  small  compared  to  the  first  order 
effect,  yielding 


=  5“  k«  4;  JJ [c„H(Q)]X.(Q)de(Q).  (3-14)’ 

With  (3-13)  we  then  obtain  at  zero  level 

•  n 

«Tti(P)  =  ^  ^  koH  TZ  ^-(P)  4;  JJ[cHH(Q)]*R„.(Q)Ar(Q) 

PM  11=1  ms— o 

and  replacing  M  by  (2-10),  the  approximated  second  order  effect  on  the 
topographic-isostatic  potential  is  simply 


«Tti(P)  =  fl<Jpolc«[cMH(P)]* 


(3-15) 


We  observe  that  under  these  approximations  the  second  order  effect  is 
independent  of  the  depth  of  compenaatioa,  ia  proportional  to  the  square  of  the 
eH*^'ivalent  rock  topog’raphy's  height,  and  ia  therefore  always  positive. 


If  we  use  the  standard  values  for  p„,  p.,  and  we  obtain  a  rule  of  thv .  b 
formula  for  the  second  order  effect  on  the  topographic-iaostatic  geoidai 
height,  using  the  Bruns  formula, 

tfNTl[e]  ^  0.2  (3-16) 

with  H*:  =  ChH,  the  equivalent  rock  topography.  Note  the  formal  similaHty  to 
the  Molodensky  effect  (HM,  p.  328) f 


Combining  first  and  second  order  effects  we  obtain 


Tti(P)  *  2e<3Dp*H*(l  +  k,  §) 


(3-17) 


from  (3-17)  we  conclude  that  the  second  order  effect,  relative  to  the  total 


potential,  increases  with  decreasing  depth  of  compensation  and  vice  versa. 
For  areas  with  extreme  topography  (both  positive  and  negative),  the  second 
order  effect  may  amount  to  about  SOX  of  the  first  order  effect  and  certainly 
cannot  be  neglected. 


4.  POWER  SPECTRUM  CONSIDERATIONS 


The  power  spectnat  (degree  variances)  of  the  topo^aphic-isostatic  potential 
in  linear  approxination  is  isplied  by  the  hamonic  coefficients  and 

of  (3-3a,b), 


Tj  :=IZ]  to  +  tIs))*.  (4-1) 

in=— fi 

Denoting  the  harmonic  coefficients  of  the  normalized  equivalent  rock 
topography  by  , 


Hj.  :=  5;  II  Ch  ^  H„.(Q)dff(Q).  (4-2) 

the  degree  variances  of  the  topographic-isostatic  potential  T^  are  related  to 
the  degree  variances  of  the  equivalent  rock  topography  Hj*, 


H*» 


n 

:=IZs;i 


(4-3) 


by 


T* 


( ^  if 

1  -  (1  -  -111 

Pn  2n  +  1  [ 

^  hJ  Jj 

i** 


(4-4) 


(cf.  Lambeckf  1979,  p.  592). 


Considering  T^  as  a  function  of  the  compensation  depth  D,  it  is  obvious  that 
the  power  spectrum  of  the  topographic-isostatic  potential  is  gaining  power 
with  increasing  compensation  depth  and  attains  a  maximum  for  D  =  R  such  that 
in  this  extreme  case  f  reduces  to  the  degree  variance  of  the  potential  of  the 
topography  for  n  >  0  (f^  =  0  in  any  case).  For  moderately  large  depths  and 
low  degrees  n  we  may  use  the  approximation  (3-5)  and  obtain  as  ratio  between 


15 


topotfraphic-isosiatic  decree  variances,  correapondintf  to  the  compensation 
depths  0|  and  D,, 

i  £i  (4-5) 

T^(Da)  Da  ' 

From  geodetic  satellite  and  surface  data  we  know  that  the  actual  anomalous 
gravitational  potential  of  the  earth  has  much  less  power  than  {T^(D  =  R)}, 
implying  that  the  rock/ocean  topography  is  isostatically  compensated  to  a 
large  extent  with  a  relatively  small  D.  A  global  standard  value  of  30  km  is 
generally  accepted  and  frequently  used  for  the  purpose  of  topographic- 
isostatic  reduction  of  surface  data. 

If  the  global  figure  for  the  compensation  depth  is  D,  we  should  expect  a  good 
agreement  between  the  power  spectrum  of  the  earth's  anootalous  gravitational 
potential  and  the  power  spectrum  Tj|(D),  provided  that  the  standard 
Airy/Heiskanen  model  describes  reality  sufficiently  well. 

This  idea  was  pursued  by  Rapp  (1982a),  using  the  simple  double  layer  model 
(linear  approximation)  and  the  DTM  data  set  consisting  of  64800  1*  *  1*  mean 
elevation  and  ocean  depth  data,  provided  by  DMAAC.  His  results  clearly 
demonstrate  that  the  standard  Airy/Heiskanen  compensation  model  with  a 
compensation  depth  of  D  =  30  km  produces  much  too  little  i>ower  over  the 
entire  frequency  range  from  n  =  2  to  180  and  is  totally  inadequate  to  explain 
the  observed  power  spectrum  of  the  anomalous  gravitational  potential.  Later 
investigations  by  Rapp  and  Rummel  (private  communication),  using  in  addition 
the  terms  of  second  and  third  order  in  H,  confirsied  the  earlier  studies. 

Guided  by  the  relation  (4-4)  or  (4-5),  considerably  larger  values  for  D  were 
suggested  by  Rapp  (1982a).  Obviously  the  best  agreesient  between  the  two 
power  spectra  was  achieved  for  a  compensation  depth  of  D  =  50  km,  yielding  a 
good  match  between  degree  50  and  ISO;  poor  results  have  been  obtained  for 
the  range  between  20  and  50  km,  but  what  is  peu^cularly  astoniahing,  also  in 
the  high  frequency  part  between  degree  150  and  180.  Moreover,  a 
compensation  depth  of  50  km  is  generally  considered  to  be  too  large  by  a 
factor  of  alsiost  2  and  is  geophysically  hard  to  Justify. 


The  South  Pacific  area  has  not  only  a  strong  appeal  to  financially  sound 
vacationers  but  also  to  hard  working  geodesists  (don’t  be  misled  by  this 
coincidence!)  because  of  extreme  gravity  field  structures.  Porsberg  (1984) 
studied  the  topographic-isostatic  geoid  in  the  Tonga  Trench  and  Tahiti  area 
using  the  recently  released  5’  «  5’  "SYNBAPS"  ocean  depth  data  set.  (I 
strongly  recommend  to  study  Forsberg's  report  -  it  is  a  beautiful  document!) 
His  results,  which  are  adso  based  on  the  standard  Airy/Heiskanen  compensation 
model,  also  suggest  much  deeper  compensation  levels  than  30  km.  Forsberg 
argues  that  at  smaller  depths  of  30  km  or  so,  conventional  isostasy  is  led  ad 
absurdum  by  deep  ocean  trenches  which  would  imply  antiroots  ending 
considerably  above  the  actual  ocean  bottom(!)  -  an  argument  in  favor  of  larger 
compensation  depths.  On  the  other  hand,  the  observed  average  thickness  of 
the  crust  below  the  ocean  bottom  is  only  about  6-8  km  -  a  fact  which  would 
favor  a  compensation  depth  smaller  than  30  km. 

Isn't  there  a  simple  way  out  of  that  isostatic  dilemma? 


4.1  Compensation  Depth  versus  Smoothing 

The  recent  studies  of  Rapp  (1982a)  and  Forsberg  (1984),  but  also  an  earlier 
work  by  Moritz  (1968)  and  Schwarz  (1976)  have  triggered  my  interest  in  this 
problem;  but  in  particular  it  was  the  simple  idea  of  Vening  Meinesz  (1939)  that 
the  actual  compensation  takes  place  regionally  rather  than  strictly  locally, 
which  made  me  to  numerically  investigate  the  problem  of  isostasy.  But  before 
we  shall  hit  the  target  in  Chapter  5,  let's  make  a  small  sidestep  to 
collocation. 

During  the  collocation  "high  noon”  period  in  the  mid-seventies,  one  of  the 
main  issues  was  how  to  deal  with  block  mean  gravity  anomalies  within  a 
homogeneous  and  isotropic  statistical  model  environment.  Numerical  integration 
of  the  covariance  function  was  prohibitive;  therefore,  an  approximation  had  to 
be  employed:  through  formally  replacing  the  block  averaging  operator  by  a 
homogeneous  and  isotropic  moving  average  operator  of  constant  weight,  the 
comfort  at  hostogeneity  and  isotropy  can  be  achieved.  According  to  the 
convolution  theorem,  the  moving  average  convolution  process  in  the  space 


domain  corresponds  to  a  simple  product  between  the  eigenvalues  of  the  moving 
average  operator  and  the  Fourier  coefficients  of  the  function  to  be  averaged 
in  the  frequency  domain.  As  far  as  collocation  is  concerned,  this  particular 
function  is  the  kernel  function  (in  least-squares  collocation  the  covariance 
function)  which,  smoothed  and  unsmoothed,  has  to  be  represented  by  a  finite 
and  as  simple  as  possible  expression  in  order  to  render  possible  a  fast 
calculation  of  linear  functionals.  Unfortunately,  the  eigenvsdues  of  the  moving 
average  operator  prevent  such  a  closed  expression  and  therefore,  they  have 
to  be  approximated  by  another  spectrum  which  behaves  properly.  Schwarz 
(1976,  p.  37  ff.)  suggested  replacing  the  moving  average  eigenvalues  with  the 
eigenvalues  of  an  upward  continuation  operator.  Applying  this  trick  means 
that  the  statistical  properties  of  mean  values  at  zero  level  are  replaced  by  the 
statistical  properties  of  point  values  at  a  certain  altitude  which  depends  on 
the  block  size  at  the  moving  average  operator.  Numerical  tests  confirmed  that 
this  approximation  wras  admissible. 

Have  we  lost  our  track?  No,  we  haven't!  The  same  idea  can  be  carried  over  to 
get,  at  least  partially,  out  of  the  iaostatic  dilemma:  assuming  that  the  Vening 
Meinesz  concept  of  regional  ieostatic  compensation  is  correct,  we  may  expect  a 
reasonably  good  solution  in  terms  of  a  power  spectrum  agreement  (between  the 
observed  anomalous  gravitational  and  the  isostasy  implied  potential)  by 
formally  replacing  the  regional  iaoeiatic  concept  of  Vening  Meineez  at 
compensation  depth  D,  by  the  strictly  local  isostatic  concept  of  Airy/Heislranen 
at  a  larger  depth  D^,  where  (Da  -  D,)  depends  on  the  characteristics  of  the 
Vening  Meinesz  smoothing. 

After  that  verbed  avalanche,  it  is  time  to  be  a  little  more  Bi)ecific:  confining 
ourselves  to  the  linear  approximation,  the  Airy/Heiskanen  power  spectrum  is 
given  by  (4-4)  which  we  repeat  here 

The  regional  Vening  Meinesz  concept  is  based  on  a  smoothing  of  the 
compensation  (  =  root/antiroot)  surface.  Since  we  don't  know  any  better,  we 
propose  democratic  smoothing  ^presented  by  a  hosiogeneous  and  isotropic 


operator  B  with  eigenvalues  p„, 


B(P,Q) 


(2n  +  l)/J„P„(co8fpq), 


(4-6a) 


1 

=  27rjB(t)P„(t)dt  (t  :=  cos  f), 

-1 


(4-6b) 


with  /Jo  =  1  to  make  the  integral  of  the  operator  B  over  the  unit  sphere  equal 
to  1  which  is  required  for  mass  balance. 

The  operator  B  will  now  be  applied  to  the  equivalent  rock  topography 
represented  by  its  harmonic  coefficients  at  the  compensation  level  only. 

Note:  the  rock/ocean  topography  remains  unchangedi  only  the  compensating 

masses  and  therefore,  the  comi>enaation  root/antiroot  surface  is  smoothed. 
This  smoothing  process  (  which  is  a  simple  convolution  in  the  space  domain)  is 
represented  by  a  multiplication  of  by  the  eigenvalue 

B  *  H*  < - ►  (4-7) 

Therefore,  the  Vening  Meiness  power  spectrum  in  its  linear  approximation  is 
given  by 


TSCD.) 


1- 

PM  2n  +  1  L 
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i*a 


(4-8) 


At  this  point  we  should  briefly  discuss  two  extreme  cases  of  sawothing: 

a)  fiff  —  1  V  n  *  Of  If  ... 

In  this  case  the  smoothing  operator  B  at  equation  (4-6a)  degenerates  into 
the  Dirac  distribution  and  will  reproduce  the  input;  consequently  there  is 
no  smoothing  involved.  Thereforef  this  operator  represents  the  standard 
Airy/Heiakanen  modeL 

b)  r  1,  =  0  V  n  =  1,  2f  ... 

In  this  case  the  smoothing  operator  B  degenerates  into  the  global  average 


op«rator  with  equal  weight;  the  compensation  masses  form  a  homogeneous 
layer  which  is  equivalent  to  a  point  mass  at  the  origin.  Therefore,  this 
operator  represents  the  extreme  case  of  the  Airy/Heiskanen  model  with 
Da  =  R. 

—  « 

It  is  evident  that  the  two  power  spectra  and  T„  coincide  if  the  condition 


(4-9) 


is  fulfilled  for  all  n.  This  implies  an  Airy/Heiskanen  compensation  depth  D,. 

Da  =  H  -  (H  -  Di)^"„  .  (4-10) 


or,  in  terms  of  radius  of  compensation,  Rf:  =  R  -  Df, 

Ra  =  RiP”  .  (4-10)' 

It  is  obvious  that  Ra  is  constant  (independent  of  the  degree  n)  if  and  only  if 
1 

p”  ~  h  =  const . , 
which  implies 

p„  =  b"  (4-11) 

with  b  <  1  because  of  the  compressing  properties  of  a  smoothing  operator. 
The  smoothing  operator  B  corresponding  to  (4-11)  is  obtained  by  an  inverse 
Founer  transform  according  to  (4-6a) 

m 

B(P,Q)  =  tIZ  (2n  +  l)b“p„(co8tp«)  (4-12a) 

^  n=o 


which  is  evidently  the  interior  Poisson  operator  (HM,  p.  35).  With  b  =  Ra/Hi 
(equation  (4-10)'  and  (4-11))  we  finally  obtain 


B(P,Q) 


R.(R?  -  RS  ) 

4wi>(P,Q) 
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(4-13) 


with  the  spatial  distance  i  defined  by 


4(P,Q)  =  (rJ  +  Ha  -  2HiHaC08fpq)’‘. 


-•  .i 

•-••A 


We  conclude:  if  the  Vening  Meineae  regional  amoothing  of  the  compensation 
surface  is  of  Poiaaon  type  at  level  Rt  with  parameter  R2,  we  can  exactly 
replace  it  by  the  standard  Airy/Heiakanen  model  with  compensation  level  Rg. 
Or  vice  veraat  usintf  a  standard  Airy/Heiakanen  model  at  depth  Da  correaponda 
exactly  to  the  use  of  a  Vening  Meineaa  model  at  depth  Di  with  the  Poiaaon 
amoothing  (4-13)  and  parameter  Da. 


Now  we  are  already  on  much  safer  ground  because  we  are  able  to  justify 
the  formal  use  of  an  unusually  large  compensation  depth.  fowever^  the 
crucial  question  is  still  oi>en:  does  the  Poisson  smoothing  reflect  physical 
reality? 


m 


5.  AN  OPTIMAL  VBNING  MEINBSZ  ISOSTATIC  MODEL 


If  we  are  talking  about  "optimal"  >  we  have  to  choose  a  norm  which  defines 
what  "good"  is.  The  choice  of  the  norm  is  not  so  difficult  if  we  are  primarily 
interested  in  the  power  spectrum  of  the  topographic-isostatic  potential.  In 
this  case  we  consider  a  solution  good  if  its  deviation  from  the  observed  power 
spectrum  of  the  anomalous  gravitational  potential  is  small.  Therefore,  it  is 
quite  natural  to  look  after  an  operator,  acting  on  the  compensation,  which 
makes  both  spectra  match  as  closely  as  possible. 

Prom  our  considerations  in  Chapter  4,  we  conclude  that  the  minimum  number 
of  parameters  must  be  2  if  we  are  after  a  compensation  smoothing  operator: 
one  peurameter  that  controls  the  level  of  compensation,  the  second  one 
controlling  the  smoothing.  (Note  that  also  the  Poisson  smoothing  operator 
(4-13)  has  two  parameters,  R|  and  Rj.) 


5.1  Estimation  of  the  parameters  of  the  saaoothing  operator 

In  linear  approximation,  the  Vening  Meinesz  topographic-isostatic  power 
spectrum  is  given  by  equation  (4-8);  the  operator’s  two  parameters  are  D  and 
b,  where  D  stands  for  depth  of  compensation  and  b  is  a  smoothing  parameter 
(note  that  depends  on  b).  As  a  matter  of  fact,  the  individual  harmonic 
coefficients  of  the  topographic-isostatic  ix>tential  are  also  controlled  in  the 
same  way, 

T„«  =  T„.(D,b).  (5-1) 

Extending  a  suggestion  of  Tscherning  (1984),  we  perform  a  least-squares 
estimation  of  the  model  p8u*ameters  such  that  the  energy  of  the  residual  field 
is  minimized. 


I  T  -  V  |*=  Bin.  (5-2) 

Since  the  low  order  harmonics  of  the  anomalous  gravitational  potential  are  to  a 
large  degree  due  to  density  disturbances  in  the  upper  mantle,  and  probably 


due  to  even  deeper  sourcesi  it  makes  no  sense  (and  would  only  disturb  the 
isostaiic  concept)  to  include  the  very  lonff  wavelength  peu't  in  the  energy 
budget. 

With  a  homogeneous  and  isotropic  smoothing  operator  B  having  eigenvalues  p„, 
the  harmonic  coefficients  of  the  topographic-isostatic  potential  are  given  by 

Tni.(D.b)  =  ^  2n  +  1  ~  ~  r1 


The  decision  of  the  model  spectrum  (j9n(^))  *>e  based  on  correlation 

■ 

considerations  between  Tjl(D,b)  and  Vj5,  requiring  that 


(5-4) 


(5-5) 


which  depend  on  the  choice  of  D.  The  figures  6.3a-d  on  page  31  show  graphs 
of  this  empirical  frequency  transfer  function  for  several  compensation 
depths  D.  As  data  the  Rapp  81  geopotential  model  and  the  worldwide  1*  *  1* 
DTM  model,  supplied  by  DMAAC,  have  been  used.  All  empirical  frequency 
transfer  functions  fi!,,  at  lecwt  for  geophysically  reasonable  compensation 
depths,  strongly  suggest  a  Gaussian  model  of  type 

^„(b)  =  e"***"’.  (5-6) 

The  solid  curves  represent  the  best  fit  of  such  a  Gaussian  model  to  the 
empirical  frequency  transfer  function. 


Having  chosen  a  model  for  it  ia  a  simple  matter  of  least-squares 

adjustment  with  2  parameters  to  solve  for  the  "best"  D  and  b.  If  we  use 
non-equal  weights  p„  which  depend  on  the  degree  like 


with  the  model  degree  variances  k„  ^nd  error  variances  of  the  harmonic 
ccefficients  of  the  anomalous  potential  modelled  by 

9%  =  ,  K  ...  constant,  (5-8) 

according  to  Jekeli  (1979,  pp.  13,  14),  we  can  even  interpret  the  least-squares 
solution  as  a  kind  of  least-squares  collocation  solution  for  the  model 
parameters  D  cmd  b  (cf.  Moritz,  1980,  pp.  160,  161).  (Here  the  data  consists  of 
the  vector  1  :=  {Vn,,}*) 


Because  of  (5-7)  and  (5-8)  the  signal  and  error  covariance  matrices  C,  and  C„ 
are  diagonal,  which  implies  that  the  energy 


2n  1 


.(D,b) 


(5-9) 


must  be  minimized  with  respect  to  D  and  b.  With  the  Taylor  linearization  of 


»D  (D=0(®),b=b(»)) 


*T,„,(D,b)  (b  -  b(«>) 

*b  (0=d(*) ,b=b(*))  (5-10) 

at  an  appropriate  Taylor  point  (for  example  D  =  30  km  and  b  according  to  the 
solid  line  in  Fig.  6.3b),  we  obtain  the  coefficients  of  the  design  matrix  A.  The 
elements  of  the  two  column  vectors  at  and  ag  consist  of  the  partiala  of  (5-10) 


'T„.(D,b) 

(o=o(**)  ,b-b(*)) 

['T„,(D,b)  I 


i 


(6-11) 


the  residual  parcuneter  vector  X  of 


XT  :=  [5D,«b].  (5-12) 

The  least-squares  solution  is  provided  by 

X  =  (aTC-‘A)~taTc-‘(1  -  !<•>)  (5-13) 

with 

C  :=  C.  +  C„  (5-14) 

and 

1<")  :=  (T„,(D<®>,b<‘*>)}  ^  (5-15) 

Finally  we  obtain  the  wanted  parameters  D  and  b  of  our  best-fittinif  Vening 
Meinesz  isostatic  model, 

D  =  +  tfD, 

b  =  +  «b. 

For  error  considerations  we  have 

B,,  =  (ATC-‘A)  ■>  (5-16) 

as  parameter  error  covariance  matrix  (for  reference  see  Moritz  (1980,  Part  B). 

5.2  Fine-tunins  of  parameter  estimation 

When  we  derived  the  least-squares  esihsation  of  the  isostatic  parameters  D 
and  b,  we  made  two  simplifications:  firstly,  if  our  Taylor  point  is  far  away 
from  reality,  we  must  iterate  the  estimation  process;  secondly,  in  our 
discussion  of  the  estinmUon  problem  and  in  all  our  derivations  we  have 

m 

feoerously  and  tacitly  assumed  that  Tn,,  of  equation  (5-3)  is  correct. 

The  first  problem  will  not  be  dealt  with  here  because  it  is  simple  and 
standard,  provided  the  Taylor  point  is  not  too  much  off  the  truth.  The 
eecond  problem  is  probably  si  ^'htly  more  delicate,  although  a  relatively  simple 
solution  is  possible: 
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Equation  (5-3)  presuppoaas  a  linaar  relation  between  the  topographic-  iaoatatic 
potential  and  the  equivalent  rock  topography.  But  according  to  (2-11)  and 
(2-12a,btc)  we  know  that  the  relation  is  highly  non-lineeu*.  Therefore,  the 
concept  of  parameter  estimation  has  to  be  modified  accordingly. 


This  modification  is  basically  due  to  second  and  higher  order  effects  which, 
we  know  from  our  discussion  in  Chapter  3,  are  small  because  the  topographic 
height  is  so  much  smaller  than  the  radius  of  the  earth.  Therefore,  an 
iteration  process  offers  itself  again. 


Let  B  be,  an  before,  the  smoothing  operator  applied  to  the  campensation 
surface,  represented  by  H* ,  yielding  the  smoothed  ccapensatioo  surface 

•I  ^  M 

H'  =  B  *  H' .  Then  H'  should  be  used  for  the  harmonic  coefficients  of  the 

•K  /  \  —  di 

ccag>en8atiQO  of  equation  (3-llb).  Note,  however,  that  H  (and  not  H)  is 

used  for  the  harsonic  coefficients  of  the  topogrwhy  because  smoothing 

is  <Hily  applied  to  the  compensating  masses.  The  problem  is  that  also  higher 
order  powers  of  H'  enter  into 


!|;(c)  ^  1 

(2n+l)(n+3) 


(l-  ^  1  5;  IJ  CHH’'»(Q)R„«(0)ds(Q) 

1^  /e 


(5-17) 


while  remains  unchanged  (equation  (3-lla)).  Since  H*  is  a  convolution  of 
B  and  H’ ,  higher  order  powers  of  H  correspond  to  convolutions  in  the  frequency 
domain.  Therefore,  the  smoothing  parameter  b  enters  in  a  rather  complicated 
manner  into  the  hi^ier  order  terms  of  T^^^*  But  fortunately,  these  terms  are 
only  small  correction  terms  to  the  leading  linear  term  (5-3);  therefore,  if  we 
solve  for  b,  it  will  be  sufficient  to  calculate  the  higher  order  terms  based  on 
the  latest  knowledge  of  b  and  to  consider  them  formally  as  constants  in  the 
subsequent  least-squares  estimation  process.  It  is  sufficient  to  iterate  this 
process  a  few  times  because  of  the  extraordinarily  rapid  convergence.  The 
entire  procedure  can  be  summarised  as  follows: 


r.  . 

r*/  s* ' 


S*  ■ 

fzl 
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hanoQic  analysis  of  low  order  powers  of  topography  according  to  (3-lla); 
correlation  study  (5-5)  yields  initial  estimates  and 

1)  L*  =  i  +  1 

smoothing  of  compensation  surface  H')  using  and 

b(‘-‘)  ; 

2)  Harmonic  analysis  of  low  order  powers  of  smoothed  compensation  surface 

according  to  (5-17)  yields  a  fiist  order  term  which  will  be 

considered  as  a  function  of  0  and  b,  and  higher  order  (correction)  terms 

m 

6T„^  which  will  be  considered  as  constants  in  the  subsequent  step; 

3)  Adjustment  (collocation)  solution  for  D  and  b  with  Taylor  point  (d(*~*), 

b^'“*^)  and  linearization  restricted  to  the  linear  term  of  step  (2); 

Stop  if  I  d(')  -  I  <  £o  and  I  b(‘^  -  b(‘"‘)  |  <  Eb,  else  go  to  (1). 

The  result  of  this  iteration  process  is  both  a  best  estioiate  for  the  depth 
of  compensation  and  for  the  parameter  of  the  most  likely  compensation 
smoothing,  and  a  set  of  harmonic  coefficients  of  the  topographic-isostatic 
potential  corresponding  to  that  Vening  Meinesz  regional  isostatic  compensation 
model. 


6.  NUMERICAL  STUDIES,  RESULTS.  CONCLUSIONS 


For  the  numerical  investitfationB  the  worldwide  digital  terrain  model  of  64800 
1*  X  1*  mean  elevations  and  ocean  bottom  depths  (but  no  information  on  ice 
coverage),  supplied  by  DMAAC  in  1979,  and,  representing  the  earth’s  gravity 
field,  the  Rapp  1981  solution  which  is  complete  to  degree  and  order  180,  have 
been  used.  Due  to  the  time  Usiitations  of  this  study  it  was  unfortunately  not 
possible  to  merge  the  recently  released  SYNBAPS  dataset,  which  consists  of 
5’  by  5’  mean  ocean  depths,  with  the  available  DTM.  Since  the  output  of  a 
system  can  hardly  be  better  than  its  input,  our  results  have  to  be  considered 
with  these  reservations  in  mind. 

As  we  mentioned  already  at  the  beginning,  the  goal  of  this  study  was  both  the 
parameter  estimation  of  the  moat  likely  iaoatatic  compensation  model  and  the 
determination  of  a  set  of  topographic-isostatic  coefficients  complete  up  to 
degree  and  order  180,  baaed  on  that  compensation  modeL  This  task  would  be 
formidable  without  support  of  the  powerful  tool  of  fast  Fourier  transform  on 
the  sphere.  Therefore,  the  unlimited  access  to  Colombo’s  (1981)  outstanding 
FFT  algorithms  "HARMIN"  for  analysis  and  "SSYMTH”  for  synthesis  was 
particularly  appreciated.  The  following  is  a  comprehensive  documentation  of 
our  numerical  studies. 

In  all  our  investigations,  the  following  parameters  have  been  used: 

density; 

Po  =  2.67  gc«-> 

p,  =  1.027  gcB-* 

Pi  =3.27  gem"* 

PN  =  5.517  gem-* 


. . .  rock  topography 
. . .  ocean  eater 
. . .  crust 


earth 


>tential  coefficient  degree  variance  aodel; 


i  i 


<n  =  (*  „i°~T.  *  =  ‘••1.  B  -  1-8 


(cf.  Rapp,  1979,  p.  10) 
potential  coefficient  error  aodel; 

<f*  =  K  =  0.0581  •  10“* 

(cf.  Jekeli,  1979,  p.  13  ff.) 

haraonic  series  coanlete  up  to  degree: 
n  =  180 

binoaial  expansion  up  to  power: 
j  =  5 

In  order  to  get  an  idea  about  the  contribution  of  first  and  higher  order  terms 
to  the  harmonic  coefficients  of  the  topographic-isostatic  potential,  we  provide 
in  Pig.  6.1  a  graph  of  the  degree  variances,  separata  for  the  linear  and 
second  order  term  as  well  as  the  composite  degree  variances.  Note  that  the 
linear  approximation  pretends  more  power  over  the  entire  frequency  rtuige 
than  the  exact  topographic-isostatic  model  actually  has.  The  contribution  of 
the  third  order  term  is  already  so  small  that  it  does  not  show  up  in  the 
plotting  window.  Therefore,  we  have  confined  our  discussion  of  higher  order 
terms  in  Chapter  3  to  the  second  order  tera  only. 


1-  .4 
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The  dependence  of  the  power  spectrum  on  the  choice  of  the  compensation 
depth  D  or,  if  Poisson  compensation  ssioothing  is  employed,  on  the 
corresponding  ssioothing  parameter  (cf.  Section  4.1),  is  illustrated  in  Pig.  6.2 
for  depths  D  =  30,  SO,  70,  and  100  km.  The  gain  of  ixiwer  with  incx*easing 
compensatk^n  depth  (or  store  compensation  ssioothing)  is  obvious. 

Frequency  transfer  functions  between  the  Rapp  81  i>»ioiiu>Wia  gravitational 
power  spectrum  and  the  topographic-laostatic  power  spectrum  in  linear 
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DECrIEf  vnHlRNCE  UOGIO  -  SCHUEl.  UNITLE5S  DECREE  VARIANCE  ILOCiO 
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Pitf.  6.1 

Degree  variances  of  the 
toi>ographic-isostatic 
poienlialt  low-order  terms 
and  total  (no  compensation 
smoothing) 
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Pig.  6.2 

Degree  variances  of  the 
topographic-isostatic 
potential  as  dependent  on 
the  depth  D  of  compensation 
(no  compensation  smoothing) 
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Pitf.  6.4a-d  Space  transfer  functions  (=  smoothintf  operators),  as  implied  by 
various  compensation  depths  D,  depending  on  spherical  distance. 


space  transfer  function  space  transfer  function 


approximation  are  given  in  Pig.  6.3a-d  for  compensation  depths  D  =  20,  30,  50, 
and  70  km.  The  graphs  strongly  suggest  a  Gaussian  model  of  the  type  (5-6). 
The  solid  lines  represent  that  beat  model  fit. 

By  an  inverse  Fourier  transform  of  the  frequency  transfer  function  p„  we 
obtain  the  corresponding  space  transfer  function  which  is  essentially  the 
compensation  smoothing  operator  B, 

B(P,Q)  =  n  (2n  +  1)^„P„(P.Q).  (6-1) 

Using  the  best  fitting  Gaussian  frequency  transfer  model,  the  smoothing 
operator  has  been  determined  for  the  most  important  part  of  its  support  and 
is  presented  in  Fig.  6.4a-d  normalized  to  B(0)  =  1.  The  scale  of  the  abscissa 
is  arc  deg.  of  spherical  distance.  These  operators  demonstrate  that  the 
smoothing  of  the  compensation  surface,  as  implied  by  our  knowledge  of  the 
global  gravitational  field  of  the  earth,  is  obviously  most  likely  of  Gaussian 
type  rather  than  of  moving  average,  box-shaped  type  or  of  Poisson  type.  The 
radius  of  smoothing  (theoretically  180  degrees,  practically  very  much  smaller) 
decreases  with  increasing  compensation  depth.  This  should  be  expected 
because  an  increase  of  compensation  depth  accounts  already  for  a  certain 
degree  of  smoothing,  although  of  Poisson  type.  For  geophysically  relevant 
compensation  depths  the  smoothing  radius  does  practically  not  exceed  2 
degrees. 

For  the  least-squares  estimation  of  the  two  model  parameters,  compensation 
depth  D  and  smoothing  parameter  b,  the  initial  (Taylor)  values 

d(«)  =  30  ka, 

(6-1) 

b^*)  =  0.0082, 

suggested  by  the  correlation  pattern,  have  been  used.  (The  smoothing 
parameter  b(^)  =  0.0082  implies  a  "correlation  length"  of  the  bell-shaped 
Gaussian  frequency  transfer  function  of  about  (degree)  n  =  100.)  Graphs  of 
both  the  frequency  and  the  corresponding  space  domain  transfer 
function(amoothing  operator)  for  these  initial  values  are  presented  in  Figs. 


SPfiCE  TRfiNSFER  FUNCTION  1=  SMOOTHING  OPERRTOR) 
BfiSED  ON  R  1  X  1  DEG.  DTM  (DMR-MODEL.  1979) 
RND  fllRT/VENlNG  HEINESZ  COMPENSATION  MODEL  WITH 
GRUSSIRN  SMOOTHING.  PARAMETER  B  =  0.0091 
DEPTH  =  24  KM.  ROCK  DENSITT  =  2.67  G/CMmm3. 
CRUST  DENSITY  =  3.27  G/CMmm3. 

WATER  DENSITT  =  1.027  G/CMx*<3 

CRLCULRTION  METHOD:  FFT  +  POWER  SERIES 


SPHERICAL  DISTANCE  (ARC. DEG.) 


Fig.  6.5  Optimal  smoothing  operator 
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SPACE  TRANSFER  FUNCTION  (UNITLESSl 


DEGREE  VRRIflNCE  (LOGIC  -  SCALE) .  UNITLESS 


DEGREE  VRRIRNCES  OF  THE  TOPOGRftPHlC-lSOSTflT IC  POTENTIRL 
BfiSED  ON  fl  I  X  1  DEG.  QTM  (DMR-HODEL,  19791 
RNO  flIRT/VENING-MElNESZ  COMPENSATION  MODEL  WITH 
ROOT/ANTIBOCTT  SMOOTHING  OF  GAUSS  TTPE,  b  -  0.0091. 
DEPTH  =  29  km.  ROCK  DENSITY  =  2.67  G/CMokS. 

*  CRUST  DENSITT  =  3.27  G/CMi<«3. 

HATER  DENSITY  -  1.027  G/CMk«3 

CRLCUlflTlON  METHOOi  FFT  ♦  POWER  SERIES 
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Fig.  6.6  Rapp  81  and  optimal  Airy/Venin*  Meineaz  degree  variances 
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The  ccaQMOisatioii  saoothing  U'  =  B  *  H*  waa  oaturally  perforaed  in  the 

m 

frequency  doaain.  The  hi|d*er  powera  of  H*  could  have  been  obtained  by  apectral 
convolution  in  the  frequency  doaain.  We  found  it  a  lot  easier  to 

m 

retransfora  the  spectrua  into  the  apace  doaain  uaing  Coloabo’s  Fourier 

synthesis  prograa  SSYNTH  and  to  calculate  all  powers  of  H*  there.  The 
binomial  exiMUision  has  been  terminated  at  power  j  =  5.  After  two  iterations  of 
the  least-squares  parameter  estimation  process  no  significant  change  of  D  and 
b  could  be  observed.  The  finally  adopted  compensation  model  parameters, 
which  yield  a  best  possible  agreement  between  gravity  field  information  and  a 
topography-isostaay  implied  model  gravity  field  on  a  global  scale,  were 


D  =  24  km 
b  =  0.0091 


(6-2) 


with  a  standard  deviation  of  less  than  2X  (!)  each.  Note  that  the  very  low 
frequency  part  of  n  <  IS  has  been  excluded  from  the  parameter  estimation 
process  in  order  to  avoid  a  strong  bias  of  the  solution  coming  from  that  part 
of  the  power  spectrum  which  has  the  least  to  do  with  isostasy.  The  "best” 
smoothing  operator  B(y)  with  b  s  0.0091,  nonaalised  to  B(0)  :  1,  is  presented 
in  Pig.  6.5  as  a  function  of  the  spherical  distance  f.  Finally  we  present  in 
Pig.  6.6  both  the  power  spectrum  of  the  Rapp  81  geopotential  solution  and  the 
power  spectrum  implied  by  the  "best"  Airy/Heiskanen-Vening  Meiness 
topographic-isoststic  model.  Because  of  the  above-oientioned  exclusion  of  the 
very  low  frequency  part  a  match  in  that  range  can  never  be  expected  (and  if 
observed,  is  purely  artificial).  The  very  long  wavelengths,  corresponding  to 
that  very  low  frequency  part,  are  to  a  great  extent  due  to  density 
distxirbances  in  the  earth's  mantle  and  are  only  sui>erimposed  by  a  relative 
little  topographic-isostatic  effect.  Therefore,  we  should  prioiarily  focus  on  the 
higher  and  highest  frequencies.  And  indeed,  in  this  range  the  agreement  of 
the  two  spectra  is  not  only  remcu’kable  -  it  is  almost  unbelievable.  At  least 
from  degree  36  on  the  siatch  is  very  good  and,  as  we  believe,  a  considerable 
improveoient  compared  with  earlier  soluttons  (compare  Rapp,  1982). 

The  compensation  ssioothing  due  to  Vening  lieiness  yields  not  only  a  better 
agreesient  with  observed  reality,  it  also  sMkes  Porsberg's  arguawnt  of 
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antiroota  with  tops  partly  above  the  ocean  bottom^  a  phenomenon  observed 
with  the  conventional  Airy/Heiskanen  model)  obsolete  because  the  compensation 
smoothing  largely  eliminates  those  singularities.  An  argument  in  favor  of  a 
smaller  compensation  depth  is  the  observed  thickness  of  the  crust  below  the 
oceans  which  is  of  the  order  of  6-8  km;  our  model  is  in  good  agreement  also 
with  this  figure.  Another  argument  which  speaks  in  favor  of  a  Vening 
Meiness  model  and  against  the  simple  Airy/Heiskanen  model  is  the  fact  that  the 
strength  of  the  earth’s  crust  is  able  to  support  a  certain  amount  of 
topographical  load  without  local  yielding,  while  in  the  1:1  model  of 
Airy/Heiskanen  free  mobility  between  vertical  mass  columns  is  presupposed  -  a 
highly  unlikely  cetse.  It  remains  an  interesting  exercise  of  applied  elasticity 
theory  to  figure  out  whether  the  Gaussian  compensation  smoothing  mth  the 
parameter  b  =  0.0091,  estimated  from  global  geodetic  evidence  only,  is 
compatible  with  the  known  elastic  parameters  of  the  crust.  Geophysicists  are 
also  invited  to  comment  on  our  geodetic  estimate  of  the  depth  of  compensation. 
In  any  case,  we  should  carefully  keep  in  mind  that  our  model  is  a  global  one 
with  only  two  parameters  -  it  can  hardly  be  simpler.  Therefore,  it  must  be 
considered  as  a  model  describing  the  average  behavior  of  the  cnist  on  a 
global  scale.  Local  deviations  from  this  global  model,  even  large  ones,  are 
possible  and  are  well-known  from  gwphysical/geological  evidence  to  exist. 
Therefore,  this  model  can  never  be  adequate  to  describe  or  even  explain 
density  patterns  which  are  due  to  plate  tectonics,  like  slabs  sliding  down  into 
the  mantle  at  plate  margins  or  the  like.  Needless  to  say,  this  was  not  the 
goal  of  our  investigation. 


Geodesists  have  been  and  still  using  the  concept  of  isostasy  for  filtering 


purposes  only,  without  bothering  very  much  about  its  geophysical  significance. 
They  are  happy  with  any  iaostatic  model  as  long  as  the  power  of  their 
measured  gravity  field  signal  is  sufficiently  diminished;  as  long  as  the 
topographic-isotatically  reduced  gravity  field  signal  is  siifflcienily  smooth, 
geodesists  can  live  with  that  model.  The  Vening  Meiness  model  proposed  here 
is  probably  not  capable  of  smoothing  gravity  field  signals  significantly  better 
than  conventional  models  which  have  been  successfully  used  in  the  past, 
although  numerical  studies  in  geodetically  pathological  areas  along  trenches 
and  continental  margins  suggest  more  optimism.  The  model  presented  here 
should  be  understood  as  a  bridge  between  geodesy  and  geophysics,  which  is 


more  stable  and  sound  than  the  ones  that  have  been  used  bjr  geodesists  in 
the  past.  It  is  a  model  which  makes  us  cheat  leas  and  therefore,  enables  us 
to  live  scientifically  more  cooifortably. 


The  harmonic  coefficients  of  the  topographic-isostatic  i>otantial  (called  TIC  85), 
which  are  based  on  the  model  parameters  (6-2)  and  derived  from  the  DMAAC 
DTM  dataset,  have  been  synthesised  using  Colombo's  Fourier  algorithm  SSYNTH, 
yielding  both  topographic-isostatic  geoidal  heights  and  gravity  anooialies  on  a 
global  1*  «  1*  grid. 

The  minimum  amd  majdmum  geoidal  height  implied  by  TIC  85  are  -10  m  and 
■»-30  m  with  a  variance  of  17  m*  which,  compared  to  the  observed  value  of 
about  900  m*,  clearly  destonstrates  the  strong  iwwer  of  non-isostatic  effects 
due  to  deep  density  sources  which  is  also  reflected  by  the  pronounced 
disagreement  between  the  TIC  85  and  the  Rapp  81  solution  in  the  very  low 
frequency  range.  It  Is  instructive  to  compare  the  minimum  and  maximum  TIC 
85  geoidal  heights  with  those  derived  from  the  rule  of  thumb  formula  (which 
does  not  presuppose  compensation  smoothing)  of  equation  (3-17):  with  the 
minimum  cmd  maximum  DTM  terrain  heights  of  -7800  m  and  ’f5900  m  equation 
(3-17)  yields  estimates  of  -12  m  and  4-25  m  which  differ  by  only  20X  from  what 
we  have  obtained.  (It  is  also  interesting  to  note  that  the  compensation 
smoothing  reduces  the  root/antiroot  maxima  by  about  30X  .)  The  histogram  of 
TIC  85  geoidal  heights  which  is  presented  in  Pig.  6.7  documents  basically  the 
distribution  of  rock  (-4)  and  ocean  (-)  masses. 

As  far  as  TIC  85  implied  gravity  anomalies  are  concerned,  we  have  minima  and 
maxima  of  -225  mgal  and  214  mgal  with  a  variance  of  516  mgal’  which  comes 
already  much  closer  to  the  observed  variance  of  1  *  k  1  *  swan  gravity 
anomalies  of  about  900  mgal*.  This  had  to  be  expected  because  gravity 
responds  much  more  to  local  density  disturbances  than  the  potential.  The 
histogram  of  TIC  85  gravity  anomalies  in  Pig.  6.8  has  a  pronounced  maximum  at 
sero  and  is  alsiost  symmetrically  bell-shaped. 

Prom  these  simple  but  instructive  statistical  data  we  can  already  anticipate  the 
Biain  features  of  the  topographic-isostatic  geoid:  it  mirrors  topography  to  a 
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great  extent  and  does  hardly  resemble  the  observed  geoid  on  a  global  scale 
because  of  the  missing  long  wavelength  information.  However,  locally  the  TIC 
geoid  comes  very  close  to  the  observed  trend-reduced  geoid.  Geoid  maps, 
produced  by  GSPP  (Sunkel,  1980)  for  various  geodetically  challenging  "hot 
spots"  on  our  earth  are  provided  as  illustration  material.  The  reader  is 
invited  to  make  a  comparison  with  the  SEASAT  -  geoid  derived  by  Rapp 
(1982b)  and  with  the  topographic-isostatic  geoid,  based  on  Airy/Heiskanen 
models  of  various  depths  and  using  the  much  more  detailed  and  accurate 
SYNBAPS  bathymetric  dataset  (S’  t  5’  compared  with  the  1*  *  1*  DMAAC  model), 
calculated  by  Porsberg  (1984)  for  some  pleasant  resort  areas  on  our  globe. 

The  TIC  8S  model  got  input  from  a  worldwide  DTM  with  a  maximum  resolution 
of  about  200  km  wavelength  and  very  poor  performance  in  some  identified 
areas  like  South  Africa,  for  example.  Because  our  operators  are  well-known  to 
be  gigo’s  (garbage  in-garbage  out),  we  strongly  suggest  to  improve  the  global 
OTM  both  with  respect  to  resolution  and  accuracy  in  order  to  make  a 
computation  of  an  even  better  model  up  to  degree  and  order  360  possible. 
Topographic-isostatic  models  of  that  high  resolution  are  not  only  important  for 
geophysical  research,  they  are  also  extremely  useful  for  the  topographic- 
isostatic  reduction  of  geodetic  data:  using  a  high  resolution  model  like  TIC  85 
or  better  as  a  global  reference,  the  entire  topographic-  isoststic  reduction 
problem  can  practically  be  reduced  to  the  processing  of  small  residuals  in 
planar  approximation.  And  for  that  purpose  we  have  again  the  very  powerful 
FFT  algorithm  at  our  disposal  which  has  been  so  successfully  applied  by 
Sideris  (1984)  for  the  evaluation  of  DTM-related  integrals,  or  alternatively,  the 
recently  published  TC-progrtuns  written  by  Forsberg(1984). 

We  consider  this  report  (which  claims  neither  to  be  complete  nor  to  be 
completely  debugged,  despite  the  heavy  use  of  a  scientific  word  processing 
system)  as  a  small  step  towards  a  better  understanding  of  our  earth’s  shell, 
which  A.  Wegener  once  compared  to  a  defendant  who  declines  to  answer.  The 
earth  scientist,  confronted  with  that  defendant,  is  the  judge  who  has  to  find 
the  truth  from  the  circumstantial  evidences. 
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